function[H1,H2]=DH(v,theta,D,k,xi,B)
%Intralyer Dirac Hamiltonian
%basis {1A,1B,2A,2B}

% k=[kx;ky];

zB=[-B(2);B(1)];

ed=pi*2.46*3.5/2.0678*1e-05;% unit is Tesla^(-1), magnetic coupling constant
 
% ed=pi*a*d/Phi, where a=2.46A, d=3.5A are lattice constants, Phi=2.0678*1e-15Wb is magnetic flux quantum

K=-4*pi/3*[1;0];

R=[cos(theta/2),-sin(theta/2);
   sin(theta/2),cos(theta/2)];

K1=R'*K;
K2=R*K;

p1=R*(k-xi*K1+ed*zB/2);
p2=R'*(k-xi*K2-ed*zB/2);

H1=[D,-v*(xi*p1(1)-1i*p1(2));
    -v*(xi*p1(1)+1i*p1(2)),-D];

H2=[D,-v*(xi*p2(1)-1i*p2(2));
    -v*(xi*p2(1)+1i*p2(2)),-D];

end